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In this work a replica exchange Monte Carlo scheme which considers an extended isobaric- 
isothermal ensemble with respect to pressure is applied to study hard spheres (HS). The idea behind 
the proposal is expanding volume instead of increasing temperature to let crowded systems char- 
acterized by dominant repulsive interactions to unblock, and so, to produce sampling from disjoint 
configurations. The method produces, in a single parallel run, the complete HS equation of state. 
Thus, the first order fluid-solid transition is captured. The obtained results well agree with previous 
calculations. This approach seems particularly useful to treat purely entropy-driven systems such 
as hard body and non-additive hard mixtures, where temperature plays a trivial role. 



I. INTRODUCTION 

The replica exchange Monte Carlo (REMC) method 
[l|, , also called parallel tempering Q , was derived to 
achieve good sampling of systems that present a free en- 
ergy landscape with many local minima [U, Q . It consists 
on simulating several replicas of the same system at dif- 
ferent thermodynamic states, and allowing for replica ex- 
changes (swap moves). Thus, it is possible to implement 
an ergodic walk through free energy barriers connecting 
disjoint configuration subspaces by defining a set of close 
enough thermodynamic states. Although it has been de- 
veloped at the end of the last century [l|, Q , its accep- 
tance is already high due to its clearness, simplicity, and 
its wide applicability. Proof of that is its employment 
to find zeolite structures [6[ , to study different conforma- 
tions of proteins 0, and to access phase equilibrium of 
many single and multicomponent systems |8l4l0l|. 

Most frequently, the REMC technique is employed to 
sample an extended canonical ensemble in temperature. 
Thus, those replicas having larger temperatures are ca- 
pable of escaping from local free energy minima, where 
the pair potential attraction of the constituting parti- 
cles plays a key role. When the free energy minima are 
mainly dictated by the entropic term, i. by the ex- 
cluded volume repulsive interactions [Ill4l3j . enlarging 
the temperature has a small effect. In other words, the 
benefits of the method become restricted. This is espe- 
cially true when dealing with hard body systems (purely 
entropy-driven systems) such as hard spheres (HS), rods, 
plates, polymers, and non-additive hard mixtures, since 
they constitute limiting cases where the pair interactions 
are repulsive only and the temperature plays a trivial 
(null) role. Thus, it is not very surprising that the REMC 
technique has not been applied yet to this kind of sys- 
tems. To do this, an alternative would be performing 
the ensemble extension in pressure instead of tempera- 
ture, to provide the particles more freedom to rearrange 
as the volume expands. This idea is tested in this work 
for HS. 

It is well known that fluid-solid transitions represent a 
challenge for computational science d [H. Most tech- 
niques which properly work for accessing liquid-gas tran- 



sitions have problems at very high densities [5|. There- 
fore, the freezing and melting points are at least difficult 
to determine [1J|. Indeed, for the HS model, simula- 
tions have recently produced an accurate determination 
of the freezin g an d melting point theoretically reported 
in the sixties [l4 [H| . That is despite the intense study 
of HS through the past decades, and the fact that the HS 
model was one of the first systems ever studied by com- 
puter simulations [l6l - [l8j . Additionally, the HS model 
shows a high density metastable branch ending at the 
random close package density which adds difficulty 
for sampling from equilibrium. 

The aim of this study is to show that the REMC can be 
successfully applied to study hard body systems. Hence, 
the REMC is used by performing a NPT ensemble exten- 
sion on pressure and applied to HS. The paper is struc- 
tured as follows. Sec. Uis this brief introduction. Sec. |TT] 
describes the employed algorithm. Results are given in 
Sec. IIIII Finally, in Sec. IIVI conclusions are drawn. 



II. NUMERICAL METHOD 

As mentioned, in the parallel tempering scheme n r 
identical replicas are considered, each following a typi- 
cal canonical simulation. However, a different temper- 
ature is set for each one of them. Thus, an extended 
ensemble can be defined so that its partition function 
is Q extended = H^QNVTi, being Q N VTi the partition 
function of ensemble i at temperature X!;, number of par- 
ticles N, and volume V. The existence of this extended 
ensemble justifies the introduction of swap trial moves 
between any two ensembles (each ensemble is sampled by 
only one replica at a time), whenever the detail balance 
condition is satisfied. If all (i, Ti)(j, Tj) (j, Ti)(i, Tj) 
swap trials have the same a priori probability of being 
performed, the swap acceptance probability becomes 

P acc =mm(l,exp[(^ - &)(£/; - [/,)]) (1) 

where Pi — l/(ksTi) is the reciprocal temperature of 
replica i, ks is the Boltzmann's constant, and Ui is the 
energy of replica i. Hence, by introducing these swap 
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trials, a particular replica seals through many tempera- 
tures allowing it to overcome free-energy barriers. Ad- 
ditionally, sampling on particular ensembles is not dis- 
turbed but enriched by the different contributions of the 
n r replicas. 

For studying systems where excluded volume interac- 
tions dominate, it may be convenient to allow the replicas 
to expand for destroying any local order. Additionally, 
in the case of a HS system (or any other purely entropy- 
driven model) an extended ensemble in temperature is 
pointless, since this variable does not affect the system 
structure. For that purpose, the extended ensemble is 
defined as Q extended = EEi Qntpu being Qntpi the 
partition function of the isobaric-isothermal ensemble of 
system i, at pressure Pi, fixed temperature T, and num- 
ber of particles N (note that the extension is in pres- 
sure; an isobaric-isotermal extension in temperature ap- 
plied to a Lennard Jones system is given by Okabe et. 
al. [Ill )■ This extended ensemble can be sampled by 
performing standard NTP simulations on each replica, 
which implies typical particle displacement trials and vol- 
ume change trials. Notwithstanding, the sampling can 
be significantly improved by introducing swap trials be- 
tween neighboring ensembles. Again, the only restriction 
is that the detail balance condition must prevail to guar- 
anty the correct sampling. One way of achieving this is 
by setting equal all a priori probabilities of choosing the 
different adjacent pairs of replicas, and accounting for 
the following acceptance probability 

P occ =mm(l,exp[/3(P - P 3 )(V - Vj)}) (2) 

where Vi — Vj is the volume difference between replicas i 
and j. It should be noted that the ensemble extension in 
pressure leads to a simple acceptance rule where energy 
terms vanish. 

For (Pi - Pj)(Vi - Vj) > 0, P acc = 1 and so, the ac- 
ceptance rule tends to order the replicas by volume size 
(lower volumes at higher pressures). For (Pj — Pj)(Vi — 
Vj) < 0, Pace depends on the absolute value of the pres- 
sure differences of the adjacent ensembles, /3\Pi — Pj\. A 
decrease of /3|P, — Pj\ leads to a larger acceptance prob- 
ability. Consequently, adjacent pressures should be close 
enough to provide large exchange acceptance rates be- 
tween neighboring ensembles. This is particularly impor- 
tant where a phase transition takes place (characterized 
by large \Vi — Vj\), which generally leads to a bottleneck 
of the swap acceptance rate. Additionally, the swap ac- 
ceptance rate also depends on the system size. Larger 
system sizes produce narrower distribution of densities 
(volumes) for a given pressure, providing smaller over- 
lap regions between adjacent ensembles. Hence, a larger 
system size leads to a decrease of the swap acceptance 
probability. Finally, in order to take a good advantage 
of the method, the replica at the lowest pressure must 
assure large jumps in configuration space, so that the 
higher pressure ensembles can be sampled from disjoin 
configurations. 

In this work, n r = 70 cubic boxes of initial side L 



are considered. These boxes are filled by randomly plac- 
ing N hard spheres of diameter a. The initial density, 
p = Nncr 3 / (6L), is set to 0.30 for all replicas. A geomet- 
rically increasing pressure, /3P, is set from approximately 
2 to 100 cr~ 3 , and arbitrarily assigned to the replicas. 
Where the fluid-solid transition is expected, intermedi- 
ate pressures are added (the total number of replicas 
equals the number of different pressures). An optimal 
allocation of replicas should lead to a constant swap ac- 
ceptance probability for all pair of adjacent ensembles 
21]. Two experiments were done, one with N — 32 and 
the other with N = 108. The simulation starts by fol- 
lowing the trial moves above described (see the appendix 
for details). 

Sampling consists on measuring densities, radial distri- 
bution functions, average number of neighbors, and the 
order parameter as a function of the pressure. The 
average number of neighbors, N n , is computed account- 
ing for all pairs having a center-center distance smaller 
than 1.2(7 (the vectors joining the centers of these pairs 
are named bonds). The order parameter Q§ is defined as 

/ rn=6 X 1 / 2 

Q°= i| E KW^)>| 2 (3) 

\ m=— 6 / 

where < Y &m (9,(f)) > is the average over all bonds and 
configurations of the spherical harmonics of the orien- 
tation angles 8 and <fi (these are the polar angles of the 
bonds measured with respect to any fixed coordinate sys- 
tem, since Qq is invariant). Qe should go to zero for a 
completely random system of a large number of points, 
following l/y/NN n /2± l/Vl3NN n [H]. 

III. RESULTS 

Figure Q] a) shows the probability density functions, 
PDFs, to find a replica at a given density for all pressures 
and for N — 32. The 70 PDFs correspond to the different 
assigned pressures. In general, the PDFs are bell-shaped 
and centered on a maximum which location depends on 
the assigned pressure. The leftmost curve corresponds 
to the lowest pressure (2.16 j3~ 1 a~ 3 ) and the rightmost 
to the highest one (100 f3~ 1 a~ 3 ). As pressure increases, 
the curves narrow and shift to the right producing larger 
densities (the narrowing is very pronounced for high pres- 
sures). The exception occurs for densities close to 0.5, 
where the PDFs split yielding bimodal distributions. At 
p ~ 0.5, the replicas produce few configurations, and 
the bimodals yield a local peak below 0.49 and another 
above 0.51. Thus, a jump on density from pf = 0.474 to 
p s = 0.520 is produced for /3P = (9.95±0.10)cr~ 3 , point- 
ing out the well known HS fluid-solid transition. The 
inset of figure [1] a) zooms in the density region around 
0.5, where the PDFs are much clearly seen. There it is 
shown the pressures that correspond to the PDFs which 
are closer to the transition. 
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FIG. 1. Probability density functions, PDFs, to find a replica 
at a given density, p. The 70 different curves correspond to the 
different assigned pressures. Fig. a) corresponds to TV =32 
and Fig. b) to TV =108. Both insets zoom in the correspond- 
ing data. 



The PDFs obtained for TV = 108 are shown in figure Q] 
b). As expected, similar trends are seen. That is, PDFs 
are bell-shaped, they narrow and shift toward larger den- 
sities for increasing pressure, and they turn bimodal for 
densities close to 0.5. Nevertheless, PDFs are higher (ap- 
proximately two times higher) and (consequently) nar- 
rower than for TV = 32. Also the bimodal distributions 
become sharper producing interpeak regions rarely vis- 
ited by the replicas. In fact, for p ~ 0.508 the PDFs 
are practically zero. In other words, the HS fluid-solid 
transition turns more evident by increasing the system 
size. As in figured] a), the inset of figure Q]b) zooms in 
the corresponding PDFs. From there it can be estimated 
the transition occurring at (3P — (10.99 ± 0.10)cr~ 3 with 
Pf = 0.487 and p s — 0.538. Thus, the transition occurs 
at a higher pressure and shifts to larger densities for in- 
creasing the system size. The gap between the fluid and 
solid densities also enlarges. 

The data obtained for small TV values can be ex- 
trapolated to estimate the HS bulk coexistence pres- 
sure, fluid density, and solid density. These are /3Pt r = 
(11.43 ±0.17)cr- 3 , p f =0.492 ± 0.004, and p s =0.545 ± 
0.004, respectively. These values are in good agreement 
with previous calculations [5j, llJ, LL5J . Figure [2] shows the 
extrapolation for the coexistence pressure, and a com- 
parison with data reported by different authors. As can 
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FIG. 2. fiP as a function of the inverse of the system size, 
1/TV. Solid symbols correspond to this work results. The 
solid symbol at 1/TV = is an extrapolation of the data and 
the solid lines are drawn to estimate the corresponding error. 
Open symbols are values reported by different authors. 



be seen, the obtained agreement is good, suggesting that 
the RECM method works properly for capturing the HS 
fluid-solid transition. 

The topmost plot of figure [3] is built by plotting the 
pressure as a function of the most frequent density for 
TV = 108. It is also shown as a red line a Pade approxi- 
mation to data obtained from the HS fluid state [23j, and 
as a blue line a fit to the HS face cubic centered (FCC) 
solid state (24|, both data series obtained by means of 
simulations. As an inset, it is shown a zoom in of the 
same data for the coexistence, where there were added 
the data obtained for TV = 32. Both curves here reported, 
for the fluid and solid states, well agree with the equa- 
tion of state given by Speedy. This confirms the good 
behavior of the REMC ensemble extension on pressure. 
Nevertheless, there is a slight deviation from the FCC 
curve of Speedy close to the transition. This may signal 
the presence of hexagonal close-packed (HCP) arrange- 
ments and even hybrid FCC-HCP structures. 

The middle and bottommost plots of figure [3] show the 
order parameter, Qq, and the number of first neighbors, 
TV„, as a function of p, respectively. The middle plot also 
shows as bullets the value of Qq for completely space- 
uncorrelated particles. As expected, Qq is small for the 
fluid region, pointing out the practical absence of angular 
order. However, it is always somewhat larger than the 
value of Qq for a random system. The difference between 
these two values diminishes for decreasing p. On the 
other hand, Qq reaches 0.5732 for f3P = 100c~ 3 , which 
is slightly lower than the Qq value of the FCC arrange- 
ment, 0.5745, and well above the corresponding value of 
the HCP structure, Q 6 = 0.4848. This signals that only 
replicas approaching the FCC lattice are allowed for the 
highest applied pressures. This is not surprising since 
108 identical spherical particles can be perfectly packed 
on a cubic box on a FCC lattice, but cannot on a HCP 
lattice. Thus, the system is being forced to promote FCC 




FIG. 3. Topmost; HS equation of state (pressure as a func- 
tion of the corresponding most frequent density) for N —108 
(o symbols). The red line corresponds to the HS fluid equa- 
tion of state of Speedy [23|] and the blue one to the HS face 
cubic centered equation of state of the same author [24|]. In- 
set; zoom in of the same plot, where □ symbols were added 
corresponding to N =32. Middle plot; Order parameter, 
as a function of p (o symbols), and the corresponding value 
for a completely random system of points 1 / \J NN„ /2 (small 
bullets). Bottommost; Number of first neighbors, N n , as a 
function of p. 



over HCP at high pressures. For lower but still over the 
coexistence pressures, Qq is close to 0.5, suggesting that 
both lattices and their hybrids contribute to the average. 
It should be noted that Qe sharply increases at the fluid- 
solid transition. Thus, it can be employed to detect any 
trace of local angular order. This was shown to be much 
more reliable than the radial distribution function peak 
that develops close to 1.5c [19]. Finally, N n monotoni- 
cally increases with p. It also shows a sharp increase at 
the coexistence, although less pronounced than for Qq. 
At large densities N n reaches 12, which is the largest 
possible HS coordination number, as it is well known. 

The radial distribution functions (RDFs) and their cor- 
responding integrals for cases a), b), c), d), and e) pointed 
out in figure [3] are plotted in figure |U Cases a) and b) 
correspond to the fluid phase, and the other three cor- 
respond to the solid phase. Case a) shows the typical 
low density liquid structure, where a relatively small con- 
tact value is developed and the second shell of neighbors 
is poorly seen. As the density increases, case b), the 
RDF shows a larger contact value (two times larger than 




FIG. 4. Radial distribution functions (black lines) and their 
integrals (red lines) for cases a), b), c), d), and e), as shown in 
figure [3] from left to right. The insets are the corresponding 
snapshots. Integrals (red lines) are scaled by a factor 1/10. 



case a)), and a well defined second shell of neighbors. 
This case corresponds to a liquid close to the coexistence. 
Slightly above the coexistence, the RDF looks like case 
c). Here a small peak appears at r/a ~ 1.5, whereas the 
valley in-between this peak and the contact one deepens. 
Other peaks also form at larger distances. For density 
values close to 0.68, case d), the RDF develops the full 
character of a crystal. That is, peaks are very high and 
narrow, and valleys turn practically zero. The integral 
(red line) of this case highlights this fact, since it shows 
a step-like behavior. The first step reaches 12, point- 
ing out the first shell coordination number (integral of 
the peak at r/a = 1), the second yields 18 (r/a = v2), 
the third 42 (r/a = y/S), and the fourth 54 (r/a = 2), 
what corresponds to the FCC structure. Nonetheless, 
a small shoulder appears at the left of the fourth peak 
(r/a = y/ll/3 ~ 1-91), suggesting the existence of few 
configurations having a HCP structure. For larger den- 
sities, case e), this shoulder disappears and a practically 
pure FCC RDF is observed. 



IV. CONCLUSIONS 

This work shows that a replica exchange Monte Carlo 
scheme can be successfully applied to study hard spheres 
at high densities. For that purpose, an extension of 
the isobaric-isothermal ensemble with respect to pres- 
sure is used. The algorithm employs standard particle 
trial displacements and volume changes together with 
replica exchanges (swap moves). These easy to imple- 
ment trials are shown to be enough for capturing the 
fluid-solid transition of hard spheres and the solid equi- 
librium branch for small systems. The obtained results 
well agree with previous calculations. The principal idea 
behind this scheme is to increase the particles mobili- 
ties by decreasing the pressure (expanding the volume), 
so that systems characterized by large excluded volume 
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contributions are able to visit disjoint configurations of 
configurational space. This approach seems particularly 
useful to deal with purely entropy-driven systems such as 
hard body and non-additive hard mixtures, where tem- 
perature plays a trivial role. 



V. APPENDIX- SIMULATION DETAILS 

Once the n r = 70 boxes are filled with the ./V spheres 
and are assigned the corresponding different pressures, 
the algorithm starts performing the trials. As mentioned, 
they are: particle displacements, volume changes, and 
swap moves. The probability for selecting a particle dis- 
placement trial (in any of the n r boxes), P^, is fixed 
to P d = n r N/(n r (N + 1) + w(n r - 1)). The probabil- 
ities for selecting a volume change trial, P v , and a swap 
trial, P s , are P v = n r /(n r (N + 1) + w{n r — 1)) and 
P s = w(n r — l)/(n r (N + 1) + w(n r — 1)). Here, w is 
a weight factor fixed to 1/50. Additionally, the proba- 
bility of performing a particle displacement trial and a 
volume change trial in a replica enlarges as it is closer 
to the fluid-solid transition pressure. These probabilities 
are 10 times larger for the central replica than for those 
having the highest and the lowest pressures. All parti- 
cles of a given replica have the same a priori probability 
of being selected to perform a displacement trial. The 
same is true for selecting a pair of adjacent replicas to 
attempt a swap move (there are n r — 1 pairs). Thus, 
a random number homogeneously distributed in [0,1] is 
generated in order to determine the type of trial to be 
performed. In case of selecting a particle displacement, 
the algorithm provides the replica and the particle for 
applying the trial. In case of a volume change trial, it 
identifies the replica; and in case of a swap trial, the al- 
gorithm gives the adjacent replicas to apply it. Next, 
another random number is generated to produce a sec- 
ond trial. If these two trials are independent the one 
another (for instance, they are particle trials on different 
replicas) the algorithm generates a third trial (note that 
these trials are not being applied yet). This procedure 
is repeated until the last trial cannot be performed inde- 
pendently of the others (for instance, a particle displace- 
ment trial on a replica in which a volume change trial 
must be previously performed) . This way, the algorithm 
have randomly selected a given number of independent 
trials to be applied on the replicas. Immediately after, 
the algorithm parallelizes (in two threads, since a dual 
core desktop is used), and the trials are done. The last 
generated trial (which was not yet performed) becomes 
now the first trial to be applied on the following series of 
trials. This procedure is followed to strictly preserve the 
detail balance condition (to build a symmetric transition 
matrix) while performing a parallelization. Verlet lists 
are employed for saving CPU time (note that the saving 
can be quite large since replicas at high pressures rarely 
update their lists). 

Sampling is not performed for the first 3.0 x 10 12 trials 




FIG. 5. a) Swap acceptance rate as a function of the ensem- 
ble pressure, b) Evolution of the densities with the number 
of trials during the initializing procedure (N = 108). As a 
reference, the data obtained from equilibrium are added as 
a solid line. The dotted lines show the random close pack- 
ing density as obtained from the square symbols and the face 
cubic centered density. 



(initializing procedure). During this process, the max- 
imum displacements of particles and maximum volume 
changes for each pressure are tuned to yield acceptance 
rates close to 0.4. Thus, particle maximum displacements 
and maximum volume changes of ensembles having high 
pressures turn smaller than those associated to ensembles 
having low pressures. Once this is done, all maximum 
displacements (one for each pressure) and maximum vol- 
ume changes (also one for each pressure) become fixed. 
1.0 x 10 13 trials are then employed to yield the data shown 
in the body of the article (equilibrium sampling). 

The acceptance rates obtained for the swap trials are 
shown in Fig. [5] a) as a function of the pressure, /3P. 
As expected, the acceptance rates for the smaller system 
(N = 32) are, in general, larger than the ones obtained 
for the larger system (N = 108). This is a consequence of 
the larger overlaps of the distributions. In both cases the 
values are above the recommended acceptance rate of 0.2 
[2lj . For j3P > PPtri the acceptance rate is practically 
constant. On the contrary, for f3P < (3Pt r , the accep- 
tance rate increases with decreasing f3P. This means that 
for the fluid region, (3P should be reduced more than geo- 
metrically for optimization purposes. For 10 < f3P < 12, 
the acceptance rate increases. This is due to the fact 
that smaller pressure differences are set between the ad- 
jacent ensembles to compensate the natural decrease of 
the acceptance rate at the fluid-solid transition. Note 
that more than the necessary replicas are added in or- 
der to decrease the error of the coexistence pressure (the 
natural decrease of the acceptance rate is overcompen- 
sated). In addition, this study is focussed on yielding a 
detailed sampling of a large /3P range. To acquire equilib- 
rium data from a high pressure system only, many fewer 
replicas would be required (the optimal swap acceptance 
rate is close to 20% when temperature is employed as the 
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thermodynamic variable of ensemble extension [21|). 

Figure [5] b) shows the evolution of the pressure versus 
density plot with the number of performed trials during 
the initializing procedure (N — 108). For ~ 10 10 tri- 
als, the sampling yields a curve at f3P > f3P tr which 
may correspond to the random close packing (RCP) 
metastable branch. There are also 5 replicas which 
reached the equilibrium state (since they crystalized, 
they were pushed towards the highest pressure region). 
Another one, producing the point laying on the dotted 
line, may correspond to a partially crystalized struc- 
ture. Assuming this well defined curve corresponds to 
the RCP branch, p rcp — 0.680 ± 0.005 is obtained from 



PP/ p 2 oc (prcp — p)^ 1 [19]. This is larger than the re- 
ported value of p rcp = 0.644 ± 0.005 [Tj3], suggesting that 
some degree of crystallization is already taking place on 
the replicas at high pressure. This is in fact confirmed by 
the Qe analysis (not shown). As the initializing process 
advances, the degree of crystallization augments and the 
high pressure curve shifts approaching the equilibrium 
branch (Qq also enlarges). For ~ 10 12 trials, the curve 
practically yields the equilibrium branch. From here on, 
only those replicas having a large degree of crystallinity 
are able to access the high pressure region. At this point, 
the initializing procedure ends and the sampling from 
equilibrium process starts. 
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